Exploración de los datos

Librerías

# Load necessary libraries
library(arm)        # Regression modeling
## Loading required package: MASS
## Loading required package: Matrix
## Loading required package: lme4
## 
## arm (Version 1.14-4, built: 2024-4-1)
## Working directory is /Users/celeste.solorzano/Desktop/Maestría Estadistica/2024/II Ciclo/Modelo Mixtos/Proyecto final
library(corrplot)   # Correlation plots
## corrplot 0.95 loaded
## 
## Attaching package: 'corrplot'
## The following object is masked from 'package:arm':
## 
##     corrplot
library(dplyr)      # Data manipulation
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)    # Data visualization
library(readxl)     # Read Excel files
library(tidyverse)  # Collection of data science tools
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ✔ readr     2.1.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ tidyr::pack()   masks Matrix::pack()
## ✖ dplyr::select() masks MASS::select()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(tidyr)      # Tidying data

Datos

# Read data
data0 <- read_excel("~/Desktop/Maestría Estadistica/2024/II Ciclo/Modelo Mixtos/Proyecto final/data.xlsx")

# Select relevant data columns
data = data0[,-c(1,3:8)]

# Convert variables to factors
data$Parte = as.factor(data$Parte)
data$Epoca = as.factor(data$Epoca)
data$Clase = as.factor(data$Clase)
# Aggregate 'Parte' into 'Seccion'
data$Seccion <- factor(
  ifelse(data$Parte %in% c("I1", "I2", "I3"), "Alta",
         ifelse(data$Parte %in% c("M1", "M2", "M3"), "Media", "Baja"))
)

# Aggregate 'Clase' into 'Conta'
data$Conta <- (
  ifelse(data$Clase %in% c("1", "2"), 0,
         ifelse(data$Clase %in% c("3", "4", "5"), 1, 0))
)

Exploración

# Descriptive analysis of the dependent variable
tabla_frecuencias <- table(data$ICA)
tabla_proporciones <- prop.table(tabla_frecuencias)
barplot(table(data$ICA), main="", xlab="Valor", ylab="Frecuencia", col="#A9C6E7")

# Correlation plot
corrplot(cor(data[,-c(1,2,28,29,30)]), method = "square", type = "lower")

# Convert ICA to factor and get the number of unique categories
data$ICA = as.factor(data$ICA)
num_categories <- length(unique(data$ICA))

# Create a color palette for the categories
formal_blue_palette <- colorRampPalette(c("#A9C6E7", "#6FA3F0", "#004B87", "#002A5C"))(num_categories)

# Function to create a bar plot using tidy evaluation
create_bar_plot <- function(x_var, x_label) {
  ggplot(data, aes(x = .data[[x_var]], fill = ICA)) +  # Use .data for tidy evaluation
    geom_bar(position = "fill", color = "black") +
    labs(title = "", x = x_label, y = "Proporción") +
    scale_y_continuous(labels = scales::percent_format()) +
    scale_fill_manual(values = formal_blue_palette) +
    theme_minimal() +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())
}

# Generate plots for each variable
plot_parte <- create_bar_plot("Parte", "Sección de la Cuenca")
plot_epoca <- create_bar_plot("Epoca", "Época")
plot_seccion <- create_bar_plot("Seccion", "Sección de la Cuenca")

# Display plots 
print(plot_parte)

print(plot_epoca)

print(plot_seccion)

# Convert response variable to numeric
data$ICA = as.numeric(data$ICA)

# Standardize numeric variables
datos_est <- data %>% mutate(across(where(is.numeric) & (3:27), scale))

# Select a sample of variables
datos_grap = datos_est[, c(1:2, 8, 22, 27:30)]

# Convert the dataframe to long format for ggplot
datos_long <- datos_grap %>%
  pivot_longer(cols = 3:4) %>%
  janitor::clean_names()  # Clean column names

# Function to create correlation plots
create_correlation_plot <- function(color_var, title) {
  ggplot(datos_long, aes(x = value, y = ica, color = !!sym(color_var))) +  
    geom_point(alpha = 0.6) +
    geom_smooth(method = "lm", se = FALSE) +
    facet_wrap(~ name) +
    labs(title = title) +
    theme_minimal()
}

# Generate correlation plots for different color variables
plot_climate <- create_correlation_plot("clase", "Correlación de Variables con Calidad de Agua por Clase")
plot_interpretation <- create_correlation_plot("interpretacion_de_calidad", "Correlación de Variables con Calidad de Agua por Interpretación de Calidad")
plot_epoca <- create_correlation_plot("epoca", "Correlación de Variables con Calidad de Agua por Época Climática")
plot_parte <- create_correlation_plot("parte", "Correlación de Variables con Calidad de Agua por Parte de la Cuenca")
plot_seccion <- create_correlation_plot("seccion", "Correlación de Variables con Calidad de Agua por Sección de la Cuenca")

# Display plots 
print(plot_climate)
## `geom_smooth()` using formula = 'y ~ x'

print(plot_interpretation)
## `geom_smooth()` using formula = 'y ~ x'

print(plot_epoca)
## `geom_smooth()` using formula = 'y ~ x'

print(plot_parte)
## `geom_smooth()` using formula = 'y ~ x'

print(plot_seccion)
## `geom_smooth()` using formula = 'y ~ x'

# Function to create heatmap
create_heatmap <- function(x_var, title) {
  ggplot(datos_long, aes_string(x = x_var, y = "epoca", fill = "ica")) +
    geom_tile() +
    scale_fill_gradient(low = "lightblue", high = "darkblue") +
    labs(title = title,
         x = "Sección de la Cuenca",
         y = "Época",
         fill = "ICA") +
    theme_minimal()
}

# Generate heatmaps for 'parte' and 'seccion'
heatmap_parte <- create_heatmap("parte", "Heatmap del Índice de Calidad de Agua por Parte")
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
heatmap_seccion <- create_heatmap("seccion", "Heatmap del Índice de Calidad de Agua por Sección")

# Display heatmaps
print(heatmap_parte)

print(heatmap_seccion)

Modelos

Librerías

# Load necessary libraries
library(MASS)
library(MCMCglmm)
## Loading required package: coda
## 
## Attaching package: 'coda'
## The following object is masked from 'package:arm':
## 
##     traceplot
## Loading required package: ape
## 
## Attaching package: 'ape'
## The following object is masked from 'package:dplyr':
## 
##     where
## The following object is masked from 'package:arm':
## 
##     balance
library(R2jags) 
## Loading required package: rjags
## Linked to JAGS 4.3.2
## Loaded modules: basemod,bugs
## 
## Attaching package: 'R2jags'
## The following object is masked from 'package:coda':
## 
##     traceplot
## The following object is masked from 'package:arm':
## 
##     traceplot
library(dplyr)
library(leaps)
library(readxl) 
library(sp)  
load.module("dic")  
## module dic loaded
load.module("glm")  
## module glm loaded

Selección de variables

# Scale selected columns (3 to 26) in the dataset
data_aux <- data %>% mutate(across(c(3:26), ~ scale(., center = TRUE, scale = FALSE)))

# Remove specific columns from the scaled dataset
data_aux1 <- data_aux[, -c(28, 29)]

# Exclude the "ICA" variable from the dataset
variables <- data_aux1[, -which(names(data_aux1) == "ICA")]

# Define the outcome variable as "ICA"
outcome <- data_aux1$ICA

# Perform best subset selection for regression
out_1 <- regsubsets(ICA ~ ., data = data_aux1)
## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax, force.in =
## force.in, : 3 linear dependencies found
## Reordering variables and trying again:
# Output the summary of the best subset selection
summary(out_1)
## Subset selection object
## Call: regsubsets.formula(ICA ~ ., data = data_aux1)
## 37 Variables  (and intercept)
##                        Forced in Forced out
## ParteF2                    FALSE      FALSE
## ParteF3                    FALSE      FALSE
## ParteI1                    FALSE      FALSE
## ParteI2                    FALSE      FALSE
## ParteI3                    FALSE      FALSE
## ParteM1                    FALSE      FALSE
## ParteM2                    FALSE      FALSE
## ParteM3                    FALSE      FALSE
## EpocaSeca                  FALSE      FALSE
## EpocaTransicion            FALSE      FALSE
## `Caudal m3/s`              FALSE      FALSE
## `T aire °C`                FALSE      FALSE
## `Vel viento m/s`           FALSE      FALSE
## `Hum Rel %`                FALSE      FALSE
## pH                         FALSE      FALSE
## `Temp agua (°C)`           FALSE      FALSE
## `OD (mg/l)`                FALSE      FALSE
## `PSO (%)`                  FALSE      FALSE
## `DBO Sol (mgO2/l)`         FALSE      FALSE
## `DBO Total (mgO2/l)`       FALSE      FALSE
## `DBO rapida`               FALSE      FALSE
## `NO3 (ug/l)`               FALSE      FALSE
## `NO2 (ug/l)`               FALSE      FALSE
## `SST (mg/l)`               FALSE      FALSE
## `Cond (µS/cm)`             FALSE      FALSE
## `NH4 (ugNH4/l)`            FALSE      FALSE
## `Alca (mg/l)`              FALSE      FALSE
## `P Tot (mgPO/l)`           FALSE      FALSE
## `P sol (mgPO/l)`           FALSE      FALSE
## `N Tot (ugN/l)`            FALSE      FALSE
## `N Org (ugN/l)`            FALSE      FALSE
## `Prom precip (mm)`         FALSE      FALSE
## `Precip acum men (mm)`     FALSE      FALSE
## Conta                      FALSE      FALSE
## `DBO Lenta`                FALSE      FALSE
## SeccionBaja                FALSE      FALSE
## SeccionMedia               FALSE      FALSE
## 1 subsets of each size up to 9
## Selection Algorithm: exhaustive
##          ParteF2 ParteF3 ParteI1 ParteI2 ParteI3 ParteM1 ParteM2 ParteM3
## 1  ( 1 ) " "     " "     " "     " "     " "     " "     " "     " "    
## 2  ( 1 ) " "     " "     " "     " "     " "     " "     " "     " "    
## 3  ( 1 ) " "     " "     " "     " "     " "     " "     " "     " "    
## 4  ( 1 ) " "     " "     " "     " "     " "     " "     " "     "*"    
## 5  ( 1 ) " "     " "     " "     " "     " "     " "     " "     " "    
## 6  ( 1 ) " "     " "     " "     " "     " "     " "     " "     "*"    
## 7  ( 1 ) " "     " "     " "     " "     " "     " "     " "     "*"    
## 8  ( 1 ) " "     " "     " "     " "     " "     " "     " "     "*"    
## 9  ( 1 ) " "     " "     " "     " "     " "     " "     "*"     "*"    
##          EpocaSeca EpocaTransicion `Caudal m3/s` `T aire °C` `Vel viento m/s`
## 1  ( 1 ) " "       " "             " "           " "         " "             
## 2  ( 1 ) " "       " "             " "           " "         " "             
## 3  ( 1 ) " "       " "             " "           " "         " "             
## 4  ( 1 ) " "       " "             " "           " "         " "             
## 5  ( 1 ) " "       " "             " "           " "         " "             
## 6  ( 1 ) " "       " "             " "           " "         " "             
## 7  ( 1 ) " "       " "             " "           " "         "*"             
## 8  ( 1 ) " "       " "             " "           " "         "*"             
## 9  ( 1 ) " "       " "             " "           "*"         "*"             
##          `Hum Rel %` pH  `Temp agua (°C)` `OD (mg/l)` `PSO (%)`
## 1  ( 1 ) " "         " " " "              " "         " "      
## 2  ( 1 ) " "         " " " "              " "         "*"      
## 3  ( 1 ) " "         " " " "              " "         "*"      
## 4  ( 1 ) " "         " " " "              " "         "*"      
## 5  ( 1 ) "*"         " " " "              " "         "*"      
## 6  ( 1 ) "*"         " " " "              " "         "*"      
## 7  ( 1 ) "*"         " " " "              " "         "*"      
## 8  ( 1 ) "*"         " " "*"              " "         "*"      
## 9  ( 1 ) "*"         " " " "              " "         "*"      
##          `DBO Sol (mgO2/l)` `DBO Total (mgO2/l)` `DBO rapida` `DBO Lenta`
## 1  ( 1 ) " "                " "                  " "          " "        
## 2  ( 1 ) " "                " "                  " "          " "        
## 3  ( 1 ) "*"                " "                  " "          " "        
## 4  ( 1 ) "*"                " "                  " "          " "        
## 5  ( 1 ) "*"                " "                  " "          " "        
## 6  ( 1 ) "*"                " "                  " "          " "        
## 7  ( 1 ) "*"                " "                  " "          " "        
## 8  ( 1 ) "*"                " "                  " "          " "        
## 9  ( 1 ) "*"                " "                  " "          " "        
##          `NO3 (ug/l)` `NO2 (ug/l)` `SST (mg/l)` `Cond (µS/cm)` `NH4 (ugNH4/l)`
## 1  ( 1 ) " "          " "          " "          " "            " "            
## 2  ( 1 ) " "          " "          " "          " "            " "            
## 3  ( 1 ) " "          " "          " "          " "            " "            
## 4  ( 1 ) " "          " "          " "          " "            " "            
## 5  ( 1 ) " "          " "          " "          " "            "*"            
## 6  ( 1 ) " "          " "          " "          " "            "*"            
## 7  ( 1 ) " "          " "          " "          " "            "*"            
## 8  ( 1 ) " "          " "          " "          " "            "*"            
## 9  ( 1 ) " "          " "          " "          " "            "*"            
##          `Alca (mg/l)` `P Tot (mgPO/l)` `P sol (mgPO/l)` `N Tot (ugN/l)`
## 1  ( 1 ) " "           " "              " "              " "            
## 2  ( 1 ) " "           " "              " "              " "            
## 3  ( 1 ) " "           " "              " "              " "            
## 4  ( 1 ) " "           " "              " "              " "            
## 5  ( 1 ) " "           " "              " "              " "            
## 6  ( 1 ) " "           " "              " "              " "            
## 7  ( 1 ) " "           " "              " "              " "            
## 8  ( 1 ) " "           " "              " "              " "            
## 9  ( 1 ) " "           " "              " "              " "            
##          `N Org (ugN/l)` `Prom precip (mm)` `Precip acum men (mm)` SeccionBaja
## 1  ( 1 ) " "             " "                " "                    " "        
## 2  ( 1 ) " "             " "                " "                    " "        
## 3  ( 1 ) " "             " "                " "                    " "        
## 4  ( 1 ) " "             " "                " "                    " "        
## 5  ( 1 ) " "             " "                " "                    " "        
## 6  ( 1 ) " "             " "                " "                    " "        
## 7  ( 1 ) " "             " "                " "                    " "        
## 8  ( 1 ) " "             " "                " "                    " "        
## 9  ( 1 ) " "             " "                " "                    " "        
##          SeccionMedia Conta
## 1  ( 1 ) " "          "*"  
## 2  ( 1 ) " "          "*"  
## 3  ( 1 ) " "          "*"  
## 4  ( 1 ) " "          "*"  
## 5  ( 1 ) " "          "*"  
## 6  ( 1 ) " "          "*"  
## 7  ( 1 ) " "          "*"  
## 8  ( 1 ) " "          "*"  
## 9  ( 1 ) " "          "*"
# Plot the results of the best subset selection
plot(out_1)

# Select variables using stepwise AIC method on a pooled linear model
modelo_seleccionado <- stepAIC(lm(ICA ~., data = data_aux1), direction = "both", trace = FALSE)

# Output the summary of the selected model
summary(modelo_seleccionado)
## 
## Call:
## lm(formula = ICA ~ Parte + `Caudal m3/s` + `Vel viento m/s` + 
##     `Hum Rel %` + pH + `OD (mg/l)` + `PSO (%)` + `DBO Total (mgO2/l)` + 
##     `DBO rapida` + `NO3 (ug/l)` + `Cond (µS/cm)` + `NH4 (ugNH4/l)` + 
##     `P Tot (mgPO/l)` + `N Tot (ugN/l)` + `N Org (ugN/l)` + `Precip acum men (mm)` + 
##     Conta, data = data_aux1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.75441 -0.34290  0.01357  0.33878  0.82511 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             4.049e+00  4.475e-01   9.048 6.07e-10 ***
## ParteF2                 5.806e-01  3.529e-01   1.645 0.110702    
## ParteF3                 4.416e-02  3.873e-01   0.114 0.910019    
## ParteI1                 1.586e-01  6.138e-01   0.258 0.797878    
## ParteI2                -1.247e-01  5.836e-01  -0.214 0.832297    
## ParteI3                 3.242e-01  5.184e-01   0.625 0.536571    
## ParteM1                 2.658e-01  4.157e-01   0.639 0.527558    
## ParteM2                -7.421e-02  4.385e-01  -0.169 0.866804    
## ParteM3                -6.608e-01  5.122e-01  -1.290 0.207210    
## `Caudal m3/s`          -4.064e-01  1.597e-01  -2.545 0.016489 *  
## `Vel viento m/s`        3.041e-01  9.254e-02   3.286 0.002663 ** 
## `Hum Rel %`             4.707e-02  1.665e-02   2.828 0.008410 ** 
## pH                      3.823e-01  2.364e-01   1.617 0.116769    
## `OD (mg/l)`            -1.774e-01  1.416e-01  -1.253 0.220226    
## `PSO (%)`              -2.081e-02  1.046e-02  -1.990 0.056044 .  
## `DBO Total (mgO2/l)`   -4.471e-02  1.137e-02  -3.934 0.000479 ***
## `DBO rapida`            7.819e-02  1.380e-02   5.665 4.01e-06 ***
## `NO3 (ug/l)`            1.174e-04  8.781e-05   1.337 0.191747    
## `Cond (µS/cm)`         -6.942e-03  2.566e-03  -2.705 0.011312 *  
## `NH4 (ugNH4/l)`         2.076e-04  4.422e-05   4.695 5.91e-05 ***
## `P Tot (mgPO/l)`        2.424e-01  1.590e-01   1.525 0.138158    
## `N Tot (ugN/l)`        -4.267e-04  2.976e-04  -1.434 0.162230    
## `N Org (ugN/l)`         4.092e-04  3.010e-04   1.360 0.184430    
## `Precip acum men (mm)`  1.071e-03  1.018e-03   1.051 0.301759    
## Conta                   2.422e+00  3.348e-01   7.234 5.77e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5459 on 29 degrees of freedom
## Multiple R-squared:  0.9694, Adjusted R-squared:  0.9442 
## F-statistic: 38.34 on 24 and 29 DF,  p-value: 4.624e-16

Modelo 1: Pooling

# Define and configure the JAGS model
jags_model_code <- "
model {
  for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau.y)           # Likelihood: y follows a normal distribution
    mu[i] <- beta0 + beta1 * x1[i] + beta2 * x2[i]  # Linear predictor
    e.y[i] <- y[i] - mu[i]                # Residuals
  }
  beta0 ~ dnorm(0, 0.001)                 # Prior for intercept
  beta1 ~ dnorm(0, 0.001)                 # Prior for slope of x1
  beta2 ~ dnorm(0, 0.001)                 # Prior for slope of x2
  tau.y <- pow(sigma.y, -2)               # Precision parameter
  sigma.y ~ dunif(0, 100)                 # Prior for standard deviation
}"
writeLines(jags_model_code, con = "modelo_jags.bug")  # Save model code to file

# Prepare data for JAGS
data_jags <- list(
  y = data$ICA,                              # Response variable
  x1 = as.numeric(data$`Temp agua (°C)`),   # Predictor variable 1
  x2 = as.numeric(data$`P sol (mgPO/l)`),    # Predictor variable 2
  n = nrow(data)                             # Number of observations
)

# Run the JAGS model
jags_model_results <- jags(
  data = data_jags,
  inits = function() list(beta0 = rnorm(1), beta1 = rnorm(1), beta2 = rnorm(1), sigma.y = runif(1, 0, 100)),  # Initial values
  parameters.to.save = c("beta0", "beta1", "beta2", "sigma.y", "e.y"),   # Parameters to monitor
  model.file = "modelo_jags.bug",   # Model file
  n.iter = 3000,                    # Total iterations
  n.burnin = 1000,                  # Burn-in period
  n.thin = 2,                       # Thinning interval
  n.chains = 3                      # Number of chains
)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 54
##    Unobserved stochastic nodes: 4
##    Total graph size: 344
## 
## Initializing model
print(jags_model_results)             # Print results
## Inference for Bugs model at "modelo_jags.bug", fit using jags,
##  3 chains, each with 3000 iterations (first 1000 discarded), n.thin = 2
##  n.sims = 3000 iterations saved. Running time = 0.176 secs
##          mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
## beta0     -8.444   2.960 -14.278 -10.368  -8.480  -6.525  -2.743 1.002  1700
## beta1      0.600   0.137   0.332   0.512   0.604   0.689   0.874 1.001  2700
## beta2      1.372   0.593   0.186   0.986   1.369   1.760   2.552 1.003   820
## e.y[1]    -1.687   0.696  -3.047  -2.141  -1.675  -1.232  -0.324 1.002  1600
## e.y[2]    -1.082   0.622  -2.285  -1.486  -1.071  -0.671   0.122 1.002  1600
## e.y[3]    -1.516   0.667  -2.830  -1.951  -1.519  -1.078  -0.210 1.002  1700
## e.y[4]    -2.416   0.486  -3.356  -2.743  -2.419  -2.094  -1.467 1.002  1700
## e.y[5]    -3.836   0.410  -4.632  -4.109  -3.841  -3.565  -3.034 1.002  1800
## e.y[6]    -5.336   0.270  -5.858  -5.513  -5.344  -5.155  -4.796 1.001  3000
## e.y[7]    -1.535   0.407  -2.357  -1.801  -1.531  -1.259  -0.737 1.001  3000
## e.y[8]     1.065   0.317   0.429   0.853   1.062   1.282   1.682 1.001  3000
## e.y[9]     1.765   0.358   1.047   1.532   1.766   2.008   2.466 1.001  3000
## e.y[10]    3.864   0.361   3.170   3.626   3.858   4.101   4.572 1.001  2100
## e.y[11]    4.044   0.389   3.296   3.783   4.036   4.302   4.805 1.001  2100
## e.y[12]    1.864   0.361   1.170   1.626   1.858   2.101   2.572 1.001  2500
## e.y[13]    2.084   0.276   1.550   1.898   2.080   2.269   2.638 1.001  3000
## e.y[14]   -0.036   0.271  -0.565  -0.217  -0.043   0.147   0.504 1.001  3000
## e.y[15]   -1.156   0.268  -1.678  -1.336  -1.164  -0.973  -0.613 1.001  3000
## e.y[16]    0.045   0.483  -0.919  -0.272   0.047   0.371   0.995 1.001  3000
## e.y[17]   -1.255   0.541  -2.366  -1.610  -1.258  -0.890  -0.197 1.001  3000
## e.y[18]   -0.435   0.577  -1.620  -0.817  -0.440  -0.044   0.690 1.001  3000
## e.y[19]    0.924   0.370   0.214   0.679   0.917   1.166   1.647 1.002  1900
## e.y[20]    0.624   0.327  -0.015   0.410   0.615   0.843   1.267 1.001  2100
## e.y[21]    0.864   0.361   0.170   0.626   0.858   1.101   1.572 1.002  1900
## e.y[22]   -0.136   0.361  -0.830  -0.374  -0.142   0.101   0.572 1.002  1900
## e.y[23]   -0.496   0.312  -1.107  -0.702  -0.506  -0.285   0.115 1.001  2200
## e.y[24]   -0.976   0.273  -1.504  -1.157  -0.980  -0.792  -0.429 1.001  3000
## e.y[25]    1.324   0.293   0.744   1.129   1.320   1.521   1.906 1.001  3000
## e.y[26]   -0.156   0.268  -0.678  -0.336  -0.164   0.027   0.387 1.001  3000
## e.y[27]   -0.216   0.268  -0.737  -0.395  -0.224  -0.032   0.323 1.001  3000
## e.y[28]   -0.903   0.401  -1.698  -1.165  -0.902  -0.633  -0.110 1.002  1400
## e.y[29]   -2.492   0.302  -3.093  -2.692  -2.492  -2.298  -1.883 1.001  2900
## e.y[30]   -0.781   0.422  -1.616  -1.058  -0.779  -0.501   0.060 1.002  1300
## e.y[31]    0.706   0.535  -0.331   0.343   0.699   1.076   1.730 1.002  1000
## e.y[32]   -0.281   0.568  -1.392  -0.664  -0.285   0.096   0.854 1.003   940
## e.y[33]    0.167   0.596  -1.035  -0.217   0.157   0.556   1.346 1.003   950
## e.y[34]   -0.609   0.460  -1.534  -0.913  -0.613  -0.311   0.309 1.002  1200
## e.y[35]   -1.345   0.441  -2.224  -1.635  -1.346  -1.055  -0.474 1.002  1200
## e.y[36]    0.109   0.575  -1.059  -0.263   0.097   0.485   1.237 1.003   980
## e.y[37]    1.667   0.607   0.453   1.263   1.668   2.073   2.839 1.003  1000
## e.y[38]    1.427   0.611   0.225   1.022   1.430   1.831   2.609 1.002  1000
## e.y[39]    2.304   0.515   1.280   1.962   2.302   2.644   3.311 1.002  1200
## e.y[40]    1.958   0.435   1.114   1.666   1.956   2.248   2.820 1.002  1200
## e.y[41]    2.143   0.308   1.531   1.939   2.137   2.349   2.750 1.001  3000
## e.y[42]    3.197   0.352   2.506   2.966   3.199   3.433   3.903 1.001  3000
## e.y[43]    3.023   0.396   2.235   2.765   3.024   3.286   3.800 1.002  1400
## e.y[44]    0.294   0.476  -0.653  -0.029   0.306   0.600   1.231 1.001  2800
## e.y[45]   -1.044   0.448  -1.935  -1.346  -1.031  -0.754  -0.168 1.001  3000
## e.y[46]   -0.510   0.567  -1.647  -0.880  -0.506  -0.130   0.612 1.002  1100
## e.y[47]   -0.666   0.558  -1.792  -1.026  -0.657  -0.298   0.452 1.003  1100
## e.y[48]   -0.680   0.545  -1.759  -1.038  -0.680  -0.321   0.378 1.002  1100
## e.y[49]   -0.671   0.322  -1.295  -0.883  -0.677  -0.451  -0.032 1.002  1900
## e.y[50]   -0.698   0.294  -1.281  -0.894  -0.699  -0.498  -0.129 1.001  3000
## e.y[51]    2.185   0.442   1.281   1.898   2.184   2.484   3.063 1.001  3000
## e.y[52]   -0.592   0.481  -1.535  -0.917  -0.576  -0.279   0.361 1.001  2500
## e.y[53]   -3.643   0.604  -4.846  -4.052  -3.635  -3.233  -2.462 1.001  3000
## e.y[54]   -0.398   0.459  -1.300  -0.707  -0.388  -0.096   0.518 1.001  2400
## sigma.y    1.924   0.196   1.589   1.785   1.908   2.040   2.363 1.001  3000
## deviance 222.200   2.935 218.548 220.085 221.490 223.655 229.525 1.001  3000
## 
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## 
## DIC info (using the rule: pV = var(deviance)/2)
## pV = 4.3 and DIC = 226.5
## DIC is an estimate of expected predictive error (lower deviance is better).
# Calculate R² for pooling
SP_pooling <- jags_model_results$BUGSoutput$sims.list   # Extract simulations from the results
r_y_pooling <- 1 - mean(apply(SP_pooling$e.y, 1, var)) / var(data$ICA)   # Calculate R² for pooling
cat(paste("R² for pooling model:", round(r_y_pooling, 4)), "\n")
## R² for pooling model: 0.3402
# Load the devtools library
library(devtools)
## Loading required package: usethis
# Install the CalvinBayes package from GitHub

devtools::install_github("rpruim/CalvinBayes")
## Skipping install of 'CalvinBayes' from a github remote, the SHA1 (dbae67a3) has not changed since last install.
##   Use `force = TRUE` to force installation
library(CalvinBayes)
## Loading required package: ggformula
## Loading required package: scales
## 
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
## The following object is masked from 'package:arm':
## 
##     rescale
## Loading required package: ggridges
## 
## New to ggformula?  Try the tutorials: 
##  learnr::run_tutorial("introduction", package = "ggformula")
##  learnr::run_tutorial("refining", package = "ggformula")
## Loading required package: bayesplot
## This is bayesplot version 1.11.1
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
##    * Does _not_ affect other ggplot2 plots
##    * See ?bayesplot_theme_set for details on theme setting
## 
## Attaching package: 'CalvinBayes'
## The following object is masked from 'package:bayesplot':
## 
##     rhat
## The following object is masked from 'package:datasets':
## 
##     HairEyeColor
# Define the parameters to analyze
parameters <- c("beta0", "beta1", "beta2")

 # Generate and display diagnostic plots 
for (param in parameters) {
  diag_mcmc(as.mcmc(jags_model_results), parName = param)
}

Modelo 2: No pooling

# Define the JAGS model with factors 'epoca' and 'parte'
jags_model_code <- "
model {
  for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau.y)                     # Likelihood: y follows a normal distribution
    mu[i] <- beta0 + beta1 * x1[i] + beta2 * x2[i] + alpha_epoca[epoca[i]] + alpha_parte[parte[i]]  # Linear predictor
    e.y[i] <- y[i] - mu[i]                          # Residuals
  }
  beta0 ~ dnorm(0, 0.001)                           # Prior for intercept
  beta1 ~ dnorm(0, 0.001)                           # Prior for slope of x1
  beta2 ~ dnorm(0, 0.001)                           # Prior for slope of x2
  
  # Priors for random effects associated with 'epoca'
  for (j in 1:J_epoca) {
    alpha_epoca[j] ~ dnorm(0, 0.001)
  }
  
  # Priors for random effects associated with 'parte'
  for (k in 1:K_parte) {
    alpha_parte[k] ~ dnorm(0, 0.001)
  }
  
  tau.y <- pow(sigma.y, -2)                         # Precision parameter
  sigma.y ~ dunif(0, 100)                           # Prior for standard deviation
}"
writeLines(jags_model_code, con = "modelo_jags.bug")  # Save model code to file

# Prepare data for JAGS
data_jags <- list(
  y = data$ICA,                                      # Response variable
  x1 = as.numeric(data$`Temp agua (°C)`),           # Predictor variable 1
  x2 = as.numeric(data$`P sol (mgPO/l)`),            # Predictor variable 2
  epoca = data$Epoca,                                # Factor for random effect 'epoca'
  parte = data$Parte,                                # Factor for random effect 'parte'
  n = nrow(data),                                    # Number of observations
  J_epoca = max(as.numeric(as.factor(data$Epoca))), # Number of levels in 'epoca'
  K_parte = max(as.numeric(as.factor(data$Parte)))   # Number of levels in 'parte'
)

# Run the JAGS model
jags_model_results <- jags(
  data = data_jags,
  inits = function() list(beta0 = rnorm(1), beta1 = rnorm(1), sigma.y = runif(1, 0, 100)),   # Initial values
  parameters.to.save = c("beta0", "beta1", "beta2", "sigma.y", "alpha_epoca", "alpha_parte", "e.y"),   # Parameters to monitor
  model.file = "modelo_jags.bug",                  # Model file path
  n.iter = 3000,                                   # Total iterations
  n.burnin = 1000,                                 # Burn-in period
  n.thin = 2,                                      # Thinning interval
  n.chains = 3                                     # Number of chains
)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 54
##    Unobserved stochastic nodes: 16
##    Total graph size: 470
## 
## Initializing model
print(jags_model_results)                           # Print results
## Inference for Bugs model at "modelo_jags.bug", fit using jags,
##  3 chains, each with 3000 iterations (first 1000 discarded), n.thin = 2
##  n.sims = 3000 iterations saved. Running time = 0.202 secs
##                mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat
## alpha_epoca[1]  -2.457  15.533 -33.474 -12.606  -2.235   7.797  27.404 1.001
## alpha_epoca[2]  -0.883  15.561 -31.928 -11.051  -0.824   9.363  28.799 1.001
## alpha_epoca[3]  -3.840  15.538 -35.005 -13.997  -3.820   6.402  25.936 1.001
## alpha_parte[1]  -0.767  10.048 -20.088  -7.515  -0.758   5.906  18.805 1.001
## alpha_parte[2]  -1.934  10.051 -21.416  -8.696  -1.937   4.620  17.586 1.001
## alpha_parte[3]  -1.228  10.034 -20.506  -7.941  -1.249   5.404  18.385 1.001
## alpha_parte[4]   0.112  10.028 -19.116  -6.568   0.070   6.781  19.895 1.001
## alpha_parte[5]  -0.179  10.038 -19.511  -6.928  -0.184   6.533  19.297 1.001
## alpha_parte[6]  -0.235  10.042 -19.327  -6.948  -0.220   6.405  19.426 1.001
## alpha_parte[7]  -1.356   9.994 -20.506  -8.043  -1.259   5.194  18.132 1.001
## alpha_parte[8]  -2.003  10.012 -21.384  -8.763  -1.946   4.549  17.670 1.001
## alpha_parte[9]  -1.731  10.050 -20.996  -8.408  -1.744   4.843  17.836 1.001
## beta0           -6.387  17.684 -41.276 -18.250  -6.339   5.333  28.065 1.001
## beta1            0.593   0.183   0.235   0.471   0.589   0.718   0.951 1.001
## beta2            3.004   0.668   1.687   2.556   2.997   3.434   4.323 1.000
## e.y[1]          -0.935   0.845  -2.608  -1.499  -0.933  -0.375   0.714 1.001
## e.y[2]          -0.077   0.839  -1.669  -0.644  -0.076   0.462   1.599 1.002
## e.y[3]          -1.115   0.843  -2.716  -1.681  -1.128  -0.556   0.553 1.001
## e.y[4]          -0.883   0.836  -2.498  -1.465  -0.878  -0.336   0.768 1.001
## e.y[5]          -1.651   0.821  -3.332  -2.194  -1.643  -1.097  -0.080 1.002
## e.y[6]          -3.405   0.758  -4.906  -3.912  -3.398  -2.896  -1.943 1.001
## e.y[7]          -0.555   0.799  -2.136  -1.094  -0.552  -0.038   1.042 1.001
## e.y[8]           3.205   0.721   1.816   2.730   3.201   3.685   4.594 1.007
## e.y[9]           3.203   0.754   1.707   2.699   3.213   3.692   4.712 1.003
## e.y[10]          0.981   0.841  -0.610   0.417   0.973   1.551   2.657 1.001
## e.y[11]          1.450   0.843  -0.249   0.895   1.452   2.003   3.086 1.002
## e.y[12]         -0.672   0.816  -2.268  -1.200  -0.689  -0.138   0.945 1.001
## e.y[13]          0.679   0.729  -0.745   0.190   0.671   1.155   2.128 1.001
## e.y[14]         -0.794   0.738  -2.278  -1.269  -0.790  -0.298   0.598 1.002
## e.y[15]         -2.184   0.786  -3.748  -2.694  -2.174  -1.671  -0.644 1.002
## e.y[16]         -1.927   0.796  -3.520  -2.462  -1.912  -1.400  -0.355 1.002
## e.y[17]         -2.056   0.753  -3.556  -2.533  -2.064  -1.568  -0.568 1.001
## e.y[18]         -1.939   0.787  -3.487  -2.450  -1.970  -1.386  -0.401 1.001
## e.y[19]         -0.386   0.821  -1.985  -0.942  -0.389   0.171   1.214 1.001
## e.y[20]         -0.392   0.837  -2.051  -0.937  -0.386   0.174   1.228 1.001
## e.y[21]         -0.099   0.781  -1.600  -0.637  -0.099   0.425   1.450 1.001
## e.y[22]          0.023   0.732  -1.412  -0.478   0.021   0.489   1.493 1.001
## e.y[23]          0.314   0.749  -1.183  -0.167   0.318   0.806   1.789 1.001
## e.y[24]         -0.432   0.771  -2.016  -0.927  -0.421   0.081   1.067 1.001
## e.y[25]          0.900   0.829  -0.723   0.348   0.893   1.457   2.503 1.002
## e.y[26]          0.593   0.785  -0.963   0.074   0.595   1.101   2.173 1.001
## e.y[27]         -0.172   0.784  -1.696  -0.700  -0.185   0.368   1.411 1.001
## e.y[28]          0.239   0.833  -1.402  -0.324   0.242   0.778   1.923 1.001
## e.y[29]         -1.469   0.818  -3.039  -1.995  -1.475  -0.926   0.181 1.001
## e.y[30]          0.782   0.819  -0.818   0.250   0.776   1.312   2.416 1.001
## e.y[31]          1.303   0.741  -0.134   0.808   1.278   1.792   2.771 1.001
## e.y[32]          0.723   0.803  -0.871   0.177   0.734   1.269   2.265 1.001
## e.y[33]          0.819   0.811  -0.881   0.302   0.830   1.351   2.393 1.001
## e.y[34]         -0.437   0.809  -2.020  -0.971  -0.443   0.095   1.144 1.001
## e.y[35]          0.020   0.756  -1.471  -0.478   0.023   0.530   1.504 1.001
## e.y[36]          0.335   0.829  -1.267  -0.199   0.352   0.886   2.003 1.001
## e.y[37]          0.533   0.810  -1.057  -0.002   0.553   1.050   2.104 1.001
## e.y[38]          0.586   0.817  -1.085   0.052   0.593   1.123   2.176 1.001
## e.y[39]          1.229   0.779  -0.319   0.721   1.236   1.729   2.726 1.001
## e.y[40]         -0.212   0.776  -1.684  -0.734  -0.213   0.275   1.354 1.001
## e.y[41]          1.273   0.743  -0.207   0.779   1.279   1.757   2.663 1.001
## e.y[42]          2.553   0.759   1.054   2.040   2.566   3.067   4.063 1.011
## e.y[43]          2.071   0.890   0.370   1.471   2.047   2.654   3.903 1.002
## e.y[44]          0.327   0.751  -1.146  -0.177   0.321   0.825   1.807 1.001
## e.y[45]         -1.975   0.756  -3.500  -2.468  -1.971  -1.487  -0.477 1.001
## e.y[46]         -0.390   0.747  -1.890  -0.865  -0.373   0.094   1.041 1.001
## e.y[47]         -0.224   0.787  -1.797  -0.746  -0.222   0.302   1.319 1.001
## e.y[48]         -0.091   0.801  -1.650  -0.628  -0.090   0.423   1.529 1.001
## e.y[49]         -0.824   0.753  -2.269  -1.352  -0.827  -0.328   0.698 1.001
## e.y[50]          0.194   0.774  -1.279  -0.335   0.186   0.728   1.712 1.001
## e.y[51]          2.630   0.820   1.020   2.111   2.651   3.163   4.206 1.001
## e.y[52]         -0.088   0.784  -1.625  -0.605  -0.101   0.430   1.476 1.001
## e.y[53]         -2.140   0.800  -3.672  -2.673  -2.154  -1.626  -0.573 1.001
## e.y[54]          0.582   0.784  -0.963   0.068   0.567   1.110   2.123 1.001
## sigma.y          1.602   0.184   1.287   1.472   1.582   1.714   2.013 1.001
## deviance       202.299   6.277 192.411 197.760 201.505 205.842 216.747 1.001
##                n.eff
## alpha_epoca[1]  3000
## alpha_epoca[2]  3000
## alpha_epoca[3]  3000
## alpha_parte[1]  3000
## alpha_parte[2]  3000
## alpha_parte[3]  3000
## alpha_parte[4]  3000
## alpha_parte[5]  3000
## alpha_parte[6]  3000
## alpha_parte[7]  3000
## alpha_parte[8]  3000
## alpha_parte[9]  3000
## beta0           3000
## beta1           3000
## beta2           3000
## e.y[1]          3000
## e.y[2]          1700
## e.y[3]          3000
## e.y[4]          3000
## e.y[5]          1900
## e.y[6]          3000
## e.y[7]          2800
## e.y[8]          3000
## e.y[9]          3000
## e.y[10]         3000
## e.y[11]         1700
## e.y[12]         2100
## e.y[13]         3000
## e.y[14]         1800
## e.y[15]         1400
## e.y[16]         1600
## e.y[17]         3000
## e.y[18]         3000
## e.y[19]         3000
## e.y[20]         3000
## e.y[21]         2200
## e.y[22]         3000
## e.y[23]         3000
## e.y[24]         2400
## e.y[25]         1600
## e.y[26]         3000
## e.y[27]         3000
## e.y[28]         3000
## e.y[29]         3000
## e.y[30]         3000
## e.y[31]         3000
## e.y[32]         3000
## e.y[33]         3000
## e.y[34]         2200
## e.y[35]         3000
## e.y[36]         3000
## e.y[37]         3000
## e.y[38]         2500
## e.y[39]         2000
## e.y[40]         3000
## e.y[41]         2700
## e.y[42]         1100
## e.y[43]         1000
## e.y[44]         3000
## e.y[45]         3000
## e.y[46]         3000
## e.y[47]         3000
## e.y[48]         2100
## e.y[49]         3000
## e.y[50]         3000
## e.y[51]         3000
## e.y[52]         2800
## e.y[53]         3000
## e.y[54]         3000
## sigma.y         2900
## deviance        2100
## 
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## 
## DIC info (using the rule: pV = var(deviance)/2)
## pV = 19.7 and DIC = 222.0
## DIC is an estimate of expected predictive error (lower deviance is better).
# Calculate R² for non-pooling model
SP_npooling <- jags_model_results$BUGSoutput$sims.list   # Extract simulations from the results
r_y_npooling <- 1 - mean(apply(SP_npooling$e.y, 1, var)) / var(data$ICA)   # Calculate R² for non-pooling model
cat(paste("R² for pooling model:", round(r_y_npooling, 4)), "\n")
## R² for pooling model: 0.5413
# Define the parameters to analyze
parameters <- c("beta0", "beta1", "beta2", "alpha_epoca[1]", "alpha_epoca[2]", "alpha_epoca[3]", "alpha_parte[1]", "alpha_parte[2]", "alpha_parte[3]", "alpha_parte[4]", "alpha_parte[5]", "alpha_parte[6]", "alpha_parte[7]", "alpha_parte[8]", "alpha_parte[9]")

 # Generate and display diagnostic plots 
for (param in parameters) {
  diag_mcmc(as.mcmc(jags_model_results), parName = param)
}

Modelo 3: Pooling parcial multinivel no anidado e intercepto variable

# Define the JAGS model
jags_model_code <- "
model {
  for (i in 1:n) {
    # Likelihood
    y[i] ~ dnorm(mu[i], tau.y)
    mu[i] <- beta0 + beta1 * x1[i] + beta2 * x2[i] + alpha_epoca[epoca[i]] + alpha_parte[parte[i]]
    
    e.y[i]  <- y[i] - mu[i]  # Residuals
  }

  # Priors for fixed effects
  beta0 ~ dnorm(0, 0.001)
  beta1 ~ dnorm(0, 0.001)
  beta2 ~ dnorm(0, 0.001)
  
  # Priors for random effects
  for (j in 1:J_epoca) {
    alpha_epoca[j] ~ dnorm(0, tau.epoca)
    e.e[j]  <- alpha_epoca[j] - beta0  # Residuals for epoch
  }
  
  for (k in 1:K_parte) {
    alpha_parte[k] ~ dnorm(0, tau.parte)
    e.p[k]  <- alpha_parte[k] - beta0   # Residuals for part
  }

  # Priors for precisions
  tau.y <- pow(sigma.y, -2)
  tau.epoca <- pow(sigma.epoca, -2)
  tau.parte <- pow(sigma.parte, -2)

  sigma.y ~ dunif(0, 100)         # Uniform prior for sigma.y
  sigma.epoca ~ dunif(0, 100)     # Uniform prior for sigma.epoca
  sigma.parte ~ dunif(0, 100)      # Uniform prior for sigma.parte
}
"

# Specify the path to save the JAGS model file
jags_model_path <- "modelo_jags.bug"

# Write the model code to a file
writeLines(jags_model_code, con = jags_model_path)

# Create a list of data for JAGS
data_jags <- list(
  y = data$ICA,                  # Dependent variable
  x1 = as.numeric(data$`Temp agua (°C)`), # Convert to numeric vector
  x2 = as.numeric(data$`P sol (mgPO/l)`),
  epoca = data$Epoca,            # Epoch index
  parte = data$Parte,            # Part index
  n = nrow(data),                # Number of observations
  J_epoca = max(as.numeric(as.factor(data$Epoca))),   # Maximum epoch index
  K_parte = max(as.numeric(as.factor(data$Parte)))    # Number of levels in part
)

# Initial values function for JAGS model parameters
inits <- function() {
  list(
    beta0 = rnorm(1),
    beta1 = rnorm(1),
    alpha_epoca = rnorm(data_jags$J_epoca),
    alpha_parte = rnorm(data_jags$K_parte),
    sigma.y = runif(1, 0, 100),
    sigma.epoca = runif(1, 0, 100),
    sigma.parte = runif(1, 0, 100)
  )
}

# Parameters to monitor during the JAGS model run
params <- c("beta0", "beta1", "beta2", "alpha_epoca", "alpha_parte", "sigma.y", "sigma.epoca", "sigma.parte", "e.e", "e.p", "e.y")

# Load and run the JAGS model
jags_model_results <- jags(
  data = data_jags,
  inits = inits,
  parameters.to.save = params,
  model.file = jags_model_path,
  n.iter = 3000,       # Total number of iterations
  n.burnin = 1000,     # Number of burn-in iterations
  n.thin = 2,          # Thinning interval (sample every second iteration)
  n.chains = 3         # Number of chains to run
)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 54
##    Unobserved stochastic nodes: 18
##    Total graph size: 486
## 
## Initializing model
# Print summary of the results from the JAGS model
print(jags_model_results)
## Inference for Bugs model at "modelo_jags.bug", fit using jags,
##  3 chains, each with 3000 iterations (first 1000 discarded), n.thin = 2
##  n.sims = 3000 iterations saved. Running time = 0.216 secs
##                mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat
## alpha_epoca[1]  -0.245   5.461 -11.679  -1.335  -0.133   0.996   9.914 1.008
## alpha_epoca[2]   1.368   5.457  -9.945   0.274   1.425   2.609  11.248 1.008
## alpha_epoca[3]  -1.433   5.465 -12.882  -2.529  -1.288  -0.210   8.564 1.008
## alpha_parte[1]   0.091   0.341  -0.570  -0.072   0.033   0.243   0.898 1.004
## alpha_parte[2]  -0.156   0.378  -1.129  -0.305  -0.066   0.035   0.463 1.001
## alpha_parte[3]   0.014   0.343  -0.715  -0.121   0.004   0.162   0.760 1.004
## alpha_parte[4]   0.165   0.383  -0.426  -0.034   0.069   0.317   1.132 1.004
## alpha_parte[5]   0.105   0.358  -0.519  -0.069   0.032   0.245   0.985 1.001
## alpha_parte[6]   0.106   0.356  -0.531  -0.060   0.035   0.239   1.031 1.002
## alpha_parte[7]  -0.045   0.338  -0.845  -0.185  -0.010   0.106   0.663 1.001
## alpha_parte[8]  -0.205   0.392  -1.256  -0.368  -0.088   0.015   0.371 1.003
## alpha_parte[9]  -0.113   0.350  -0.983  -0.254  -0.035   0.059   0.496 1.002
## beta0           -5.879   6.142 -17.104  -8.454  -5.957  -3.318   6.037 1.007
## beta1            0.447   0.132   0.193   0.359   0.444   0.531   0.710 1.001
## beta2            2.304   0.577   1.211   1.908   2.291   2.670   3.480 1.001
## e.e[1]           5.634  11.261 -17.174   2.414   5.835   9.102  25.649 1.008
## e.e[2]           7.247  11.235 -15.456   4.017   7.439  10.683  27.616 1.008
## e.e[3]           4.446  11.278 -18.915   1.202   4.736   7.959  24.228 1.008
## e.p[1]           5.970   6.138  -5.847   3.371   6.055   8.560  17.147 1.007
## e.p[2]           5.723   6.107  -6.344   3.240   5.758   8.286  16.745 1.007
## e.p[3]           5.893   6.127  -5.989   3.323   6.003   8.513  16.825 1.007
## e.p[4]           6.044   6.202  -5.996   3.421   6.084   8.639  17.437 1.007
## e.p[5]           5.984   6.203  -6.049   3.287   6.024   8.584  17.369 1.007
## e.p[6]           5.985   6.206  -6.084   3.314   6.068   8.598  17.307 1.007
## e.p[7]           5.834   6.166  -6.173   3.266   5.884   8.440  16.995 1.007
## e.p[8]           5.674   6.150  -6.435   3.111   5.743   8.321  16.865 1.007
## e.p[9]           5.766   6.137  -6.208   3.220   5.778   8.354  16.800 1.007
## e.y[1]          -0.969   0.672  -2.337  -1.410  -0.960  -0.521   0.319 1.001
## e.y[2]          -0.236   0.625  -1.462  -0.656  -0.231   0.176   0.988 1.001
## e.y[3]          -1.078   0.634  -2.329  -1.474  -1.072  -0.659   0.097 1.001
## e.y[4]          -1.598   0.554  -2.687  -1.955  -1.600  -1.255  -0.497 1.001
## e.y[5]          -2.751   0.563  -3.770  -3.133  -2.785  -2.401  -1.527 1.002
## e.y[6]          -3.961   0.490  -4.940  -4.275  -3.975  -3.653  -2.949 1.002
## e.y[7]          -0.060   0.600  -1.302  -0.438  -0.040   0.329   1.072 1.001
## e.y[8]           2.634   0.528   1.599   2.289   2.621   2.971   3.736 1.002
## e.y[9]           3.241   0.544   2.127   2.891   3.249   3.597   4.301 1.001
## e.y[10]          1.854   0.596   0.602   1.459   1.879   2.264   2.966 1.001
## e.y[11]          2.049   0.593   0.858   1.667   2.059   2.449   3.202 1.001
## e.y[12]         -0.086   0.583  -1.279  -0.467  -0.069   0.301   1.046 1.001
## e.y[13]          0.483   0.504  -0.481   0.146   0.486   0.805   1.476 1.001
## e.y[14]         -1.446   0.527  -2.435  -1.802  -1.457  -1.113  -0.329 1.003
## e.y[15]         -2.628   0.512  -3.591  -2.977  -2.631  -2.309  -1.583 1.001
## e.y[16]         -1.174   0.589  -2.371  -1.558  -1.141  -0.786  -0.044 1.001
## e.y[17]         -2.151   0.595  -3.275  -2.548  -2.156  -1.780  -0.961 1.001
## e.y[18]         -1.454   0.606  -2.662  -1.849  -1.450  -1.061  -0.262 1.001
## e.y[19]          0.512   0.564  -0.680   0.165   0.528   0.894   1.583 1.002
## e.y[20]          0.349   0.538  -0.778   0.009   0.367   0.720   1.361 1.001
## e.y[21]          0.527   0.542  -0.617   0.175   0.546   0.890   1.547 1.002
## e.y[22]         -0.322   0.518  -1.335  -0.666  -0.321   0.021   0.703 1.001
## e.y[23]         -0.431   0.528  -1.390  -0.793  -0.447  -0.094   0.697 1.001
## e.y[24]         -0.881   0.496  -1.831  -1.217  -0.883  -0.570   0.123 1.001
## e.y[25]          1.139   0.511   0.062   0.816   1.147   1.493   2.116 1.001
## e.y[26]          0.028   0.504  -0.949  -0.313   0.016   0.347   1.058 1.001
## e.y[27]         -0.187   0.492  -1.176  -0.502  -0.179   0.139   0.743 1.002
## e.y[28]          0.662   0.600  -0.535   0.267   0.680   1.068   1.817 1.001
## e.y[29]         -1.042   0.543  -2.141  -1.380  -1.026  -0.685  -0.002 1.001
## e.y[30]          0.871   0.613  -0.326   0.465   0.866   1.279   2.071 1.001
## e.y[31]          1.257   0.531   0.194   0.911   1.267   1.614   2.282 1.001
## e.y[32]          0.486   0.572  -0.645   0.102   0.490   0.866   1.606 1.002
## e.y[33]          0.915   0.596  -0.290   0.510   0.924   1.309   2.076 1.001
## e.y[34]          0.256   0.579  -0.931  -0.102   0.275   0.653   1.332 1.001
## e.y[35]         -0.277   0.528  -1.286  -0.630  -0.280   0.060   0.798 1.002
## e.y[36]          0.803   0.615  -0.450   0.413   0.817   1.222   1.964 1.001
## e.y[37]          0.819   0.620  -0.397   0.396   0.828   1.238   2.011 1.001
## e.y[38]          0.701   0.628  -0.494   0.273   0.697   1.118   1.944 1.001
## e.y[39]          1.381   0.573   0.249   0.995   1.381   1.755   2.498 1.001
## e.y[40]          0.039   0.597  -1.162  -0.348   0.043   0.427   1.184 1.001
## e.y[41]          0.847   0.517  -0.119   0.501   0.835   1.178   1.922 1.002
## e.y[42]          2.182   0.500   1.241   1.848   2.167   2.500   3.198 1.002
## e.y[43]          1.867   0.553   0.812   1.508   1.854   2.227   2.978 1.001
## e.y[44]         -0.333   0.573  -1.373  -0.720  -0.354   0.023   0.862 1.002
## e.y[45]         -1.959   0.518  -2.974  -2.310  -1.956  -1.623  -0.949 1.002
## e.y[46]         -0.195   0.579  -1.322  -0.587  -0.184   0.183   0.953 1.001
## e.y[47]         -0.228   0.587  -1.410  -0.627  -0.226   0.168   0.948 1.001
## e.y[48]          0.077   0.592  -1.096  -0.321   0.088   0.468   1.257 1.001
## e.y[49]         -0.723   0.510  -1.726  -1.050  -0.720  -0.378   0.245 1.001
## e.y[50]         -0.274   0.516  -1.228  -0.619  -0.300   0.048   0.834 1.001
## e.y[51]          2.655   0.571   1.523   2.284   2.643   3.021   3.781 1.001
## e.y[52]          0.169   0.564  -0.961  -0.207   0.180   0.530   1.308 1.001
## e.y[53]         -2.508   0.663  -3.766  -2.945  -2.519  -2.084  -1.208 1.001
## e.y[54]          0.403   0.553  -0.683   0.038   0.407   0.758   1.498 1.001
## sigma.epoca      6.501  10.346   0.779   1.694   2.938   6.174  38.832 1.001
## sigma.parte      0.350   0.287   0.004   0.134   0.288   0.497   1.070 1.161
## sigma.y          1.570   0.166   1.292   1.451   1.556   1.676   1.927 1.001
## deviance       200.126   4.394 192.984 197.144 199.519 202.534 210.752 1.002
##                n.eff
## alpha_epoca[1]  3000
## alpha_epoca[2]  3000
## alpha_epoca[3]  3000
## alpha_parte[1]  1800
## alpha_parte[2]  3000
## alpha_parte[3]  1100
## alpha_parte[4]   780
## alpha_parte[5]  3000
## alpha_parte[6]  1100
## alpha_parte[7]  3000
## alpha_parte[8]  1500
## alpha_parte[9]  3000
## beta0           3000
## beta1           3000
## beta2           3000
## e.e[1]          3000
## e.e[2]          3000
## e.e[3]          3000
## e.p[1]          2700
## e.p[2]          3000
## e.p[3]          2700
## e.p[4]          2800
## e.p[5]          2900
## e.p[6]          3000
## e.p[7]          3000
## e.p[8]          3000
## e.p[9]          2700
## e.y[1]          3000
## e.y[2]          3000
## e.y[3]          3000
## e.y[4]          3000
## e.y[5]          1000
## e.y[6]          1800
## e.y[7]          2900
## e.y[8]          1200
## e.y[9]          2800
## e.y[10]         3000
## e.y[11]         3000
## e.y[12]         3000
## e.y[13]         3000
## e.y[14]          930
## e.y[15]         3000
## e.y[16]         3000
## e.y[17]         3000
## e.y[18]         3000
## e.y[19]         1800
## e.y[20]         3000
## e.y[21]         2200
## e.y[22]         3000
## e.y[23]         2200
## e.y[24]         3000
## e.y[25]         3000
## e.y[26]         3000
## e.y[27]         2300
## e.y[28]         3000
## e.y[29]         3000
## e.y[30]         3000
## e.y[31]         3000
## e.y[32]         1700
## e.y[33]         3000
## e.y[34]         3000
## e.y[35]         1900
## e.y[36]         3000
## e.y[37]         3000
## e.y[38]         3000
## e.y[39]         3000
## e.y[40]         3000
## e.y[41]         1100
## e.y[42]         3000
## e.y[43]         3000
## e.y[44]         1700
## e.y[45]         3000
## e.y[46]         3000
## e.y[47]         3000
## e.y[48]         3000
## e.y[49]         3000
## e.y[50]         2900
## e.y[51]         3000
## e.y[52]         3000
## e.y[53]         3000
## e.y[54]         3000
## sigma.epoca     3000
## sigma.parte       37
## sigma.y         3000
## deviance        1100
## 
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## 
## DIC info (using the rule: pV = var(deviance)/2)
## pV = 9.6 and DIC = 209.8
## DIC is an estimate of expected predictive error (lower deviance is better).
# Calculate R²
SP_model3 <- jags_model_results$BUGSoutput$sims.list   # Extract simulations from the results
r_y_model3 <- 1 - mean(apply(SP_model3$e.y, 1, var)) / var(data$ICA)   # Calculate R² 
cat(paste("R² for model 3:", round(r_y_model3, 4)), "\n")
## R² for model 3: 0.561
# Extract posterior samples and calculate omega
posterior_samples <- jags_model_results$BUGSoutput$sims.list
omega_model3 <- pmin((apply(posterior_samples$e.y, 2, sd) / mean(posterior_samples$sigma.y))^2, 1)
cat("Omega:\n")
## Omega:
plot(
  omega_model3,
  type = "l", col = "blue", lwd = 2,
  main = "", xlab = "Index", ylab = "Omega"
)

# Define the parameters to analyze
parameters <- c("beta0", "beta1", "beta2", "alpha_epoca[1]", "alpha_epoca[2]", "alpha_epoca[3]", "alpha_parte[1]", "alpha_parte[2]", "alpha_parte[3]", "alpha_parte[4]", "alpha_parte[5]", "alpha_parte[6]", "alpha_parte[7]", "alpha_parte[8]", "alpha_parte[9]")

 # Generate and display diagnostic plots 
for (param in parameters) {
  diag_mcmc(as.mcmc(jags_model_results), parName = param)
}

Modelo 4: Pooling parcial multinivel con efectos aleatorios a nivel de pendientes

# Define the JAGS model
jags_model_code <- "
model {
  for (i in 1:n) {
    # Likelihood
    y[i] ~ dnorm(mu[i], tau.y)
    mu[i] <- beta0 + beta1[epoca[i]] * x1[i] + beta2[epoca[i]] * x2[i]
    e.y[i]  <- y[i] - mu[i]  # Residuals
  }

  # Priors for fixed effects
  beta0 ~ dnorm(0, 0.001)
  
  # Priors for precision of the likelihood
  tau.y <- pow(sigma.y, -2)
  sigma.y ~ dunif(0, 100)
  
  # Priors for random effects
  for (j in 1:J_epoca) {
    beta1[j] ~ dnorm(m1, t1)
    beta2[j] ~ dnorm(m2, t2)
  }
  
  # Hyperparameters for random effects
  m1 ~ dnorm(0, 0.0001)
  m2 ~ dnorm(0, 0.0001)
  
  t1 <- pow(s1, -2)
  s1 ~ dunif(0, 100)
  
  t2 <- pow(s2, -2)
  s2 ~ dunif(0, 100)
}
"

# Specify the path to save the JAGS model file
jags_model_path <- "modelo_jags.bug"

# Write the model code to a file
writeLines(jags_model_code, con = jags_model_path)

# Create a list of data for JAGS
data_jags <- list(
  y = data$ICA,                  # Dependent variable
  x1 = as.numeric(data$`Temp agua (°C)`), # Independent variable 1
  x2 = as.numeric(data$`P sol (mgPO/l)`), # Independent variable 2
  epoca = data$Epoca,            # Epoch index
  n = nrow(data),                # Number of observations
  J_epoca = max(as.numeric(as.factor(data$Epoca))) # Maximum epoch index
)

# Initial values function for JAGS model parameters
inits <- function() {
  list(
    sigma.y = runif(1, 0, 100),   # Initial value for sigma.y
    s1 = runif(1, 0, 100),        # Initial value for s1
    s2 = runif(1, 0, 100)         # Initial value for s2
  )
}

# Parameters to monitor during the JAGS model run
params <- c("beta0", "beta1", "beta2", "sigma.y", "s1", "s2", "m1", "m2", "e.y")

# Load and run the JAGS model
jags_model_results <- jags(
  data = data_jags,
  inits = inits,
  parameters.to.save = params,
  model.file = jags_model_path,
  n.iter = 10000,       # Total number of iterations
  n.burnin = 3000,     # Number of burn-in iterations
  n.thin = 4,          # Thinning interval (sample every second iteration)
  n.chains = 3         # Number of chains to run
)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 54
##    Unobserved stochastic nodes: 12
##    Total graph size: 424
## 
## Initializing model
# Print summary of the results from the JAGS model
print(jags_model_results)
## Inference for Bugs model at "modelo_jags.bug", fit using jags,
##  3 chains, each with 10000 iterations (first 3000 discarded), n.thin = 4
##  n.sims = 5250 iterations saved. Running time = 0.342 secs
##          mu.vect sd.vect    2.5%     25%     50%     75%   97.5%  Rhat n.eff
## beta0     -6.031   2.382 -10.305  -7.858  -6.093  -4.344  -1.317 1.027    79
## beta1[1]   0.439   0.114   0.212   0.357   0.442   0.523   0.652 1.025    93
## beta1[2]   0.526   0.112   0.307   0.445   0.529   0.610   0.730 1.026    84
## beta1[3]   0.399   0.124   0.158   0.308   0.404   0.492   0.624 1.023    90
## beta2[1]   2.405   0.856   0.843   1.841   2.352   2.926   4.245 1.001  5200
## beta2[2]   1.895   0.809   0.265   1.375   1.923   2.417   3.468 1.001  2800
## beta2[3]   2.156   0.771   0.656   1.654   2.147   2.654   3.727 1.001  3300
## e.y[1]    -1.164   0.591  -2.335  -1.570  -1.160  -0.756  -0.029 1.017   130
## e.y[2]    -0.459   0.549  -1.549  -0.830  -0.451  -0.085   0.598 1.014   150
## e.y[3]    -1.250   0.556  -2.326  -1.631  -1.252  -0.864  -0.188 1.019   110
## e.y[4]    -1.849   0.440  -2.702  -2.141  -1.852  -1.555  -0.990 1.012   180
## e.y[5]    -3.128   0.405  -3.918  -3.399  -3.128  -2.860  -2.324 1.007   300
## e.y[6]    -4.125   0.425  -4.959  -4.412  -4.128  -3.840  -3.271 1.002  1600
## e.y[7]     0.077   0.577  -1.042  -0.327   0.077   0.466   1.229 1.009   250
## e.y[8]     2.476   0.492   1.517   2.141   2.475   2.806   3.451 1.005   430
## e.y[9]     3.277   0.533   2.235   2.909   3.278   3.637   4.337 1.007   310
## e.y[10]    2.380   0.504   1.403   2.041   2.373   2.716   3.373 1.010   230
## e.y[11]    2.538   0.523   1.533   2.187   2.532   2.892   3.569 1.011   200
## e.y[12]    0.380   0.504  -0.597   0.041   0.373   0.716   1.373 1.010   210
## e.y[13]    0.696   0.443  -0.155   0.398   0.691   0.993   1.569 1.004   570
## e.y[14]   -1.409   0.438  -2.250  -1.701  -1.414  -1.116  -0.553 1.003   740
## e.y[15]   -2.514   0.433  -3.342  -2.804  -2.519  -2.224  -1.669 1.003  1000
## e.y[16]   -1.093   0.495  -2.044  -1.426  -1.099  -0.761  -0.106 1.005   430
## e.y[17]   -2.356   0.526  -3.346  -2.709  -2.368  -1.997  -1.316 1.007   310
## e.y[18]   -1.514   0.546  -2.547  -1.885  -1.525  -1.144  -0.424 1.008   260
## e.y[19]    0.634   0.470  -0.311   0.330   0.642   0.949   1.563 1.008   270
## e.y[20]    0.414   0.451  -0.480   0.123   0.418   0.716   1.291 1.006   400
## e.y[21]    0.590   0.466  -0.351   0.289   0.598   0.904   1.508 1.008   290
## e.y[22]   -0.410   0.466  -1.351  -0.711  -0.402  -0.096   0.508 1.008   290
## e.y[23]   -0.674   0.444  -1.562  -0.965  -0.670  -0.374   0.189 1.005   480
## e.y[24]   -1.025   0.432  -1.889  -1.307  -1.023  -0.731  -0.199 1.002  1500
## e.y[25]    1.194   0.437   0.318   0.905   1.199   1.488   2.041 1.004   690
## e.y[26]   -0.157   0.432  -1.016  -0.440  -0.157   0.138   0.664 1.001  2900
## e.y[27]   -0.201   0.432  -1.067  -0.487  -0.202   0.096   0.618 1.001  3700
## e.y[28]    0.637   0.679  -0.654   0.188   0.636   1.088   1.987 1.001  5200
## e.y[29]   -1.070   0.551  -2.135  -1.447  -1.076  -0.696   0.026 1.001  4700
## e.y[30]    0.774   0.703  -0.572   0.311   0.768   1.239   2.173 1.001  5200
## e.y[31]    1.179   0.568   0.024   0.807   1.185   1.545   2.286 1.003  1100
## e.y[32]    0.334   0.620  -0.906  -0.059   0.337   0.733   1.538 1.001  5200
## e.y[33]    0.901   0.644  -0.374   0.483   0.904   1.318   2.154 1.002  2400
## e.y[34]    0.417   0.515  -0.589   0.074   0.420   0.749   1.435 1.003   780
## e.y[35]   -0.385   0.497  -1.365  -0.712  -0.382  -0.064   0.607 1.002  1300
## e.y[36]    0.919   0.621  -0.307   0.516   0.917   1.323   2.128 1.002  1600
## e.y[37]    0.823   0.748  -0.667   0.328   0.835   1.316   2.268 1.001  5200
## e.y[38]    0.612   0.757  -0.903   0.119   0.625   1.111   2.065 1.001  5200
## e.y[39]    1.382   0.634   0.129   0.962   1.389   1.797   2.617 1.001  5200
## e.y[40]    0.382   0.688  -0.977  -0.058   0.360   0.830   1.749 1.002  1900
## e.y[41]    0.820   0.453  -0.058   0.517   0.809   1.126   1.698 1.001  5200
## e.y[42]    2.076   0.411   1.275   1.798   2.072   2.342   2.889 1.002  2000
## e.y[43]    1.962   0.489   1.004   1.630   1.965   2.284   2.944 1.002  2300
## e.y[44]   -0.640   0.518  -1.636  -0.992  -0.641  -0.301   0.394 1.004   640
## e.y[45]   -2.046   0.472  -2.943  -2.361  -2.054  -1.740  -1.106 1.004   530
## e.y[46]    0.011   0.651  -1.231  -0.420   0.000   0.425   1.346 1.006   420
## e.y[47]   -0.078   0.661  -1.334  -0.516  -0.096   0.340   1.281 1.004   550
## e.y[48]    0.247   0.719  -1.093  -0.244   0.229   0.714   1.744 1.001  5200
## e.y[49]   -0.815   0.530  -1.871  -1.158  -0.807  -0.452   0.198 1.001  5200
## e.y[50]   -0.496   0.440  -1.359  -0.788  -0.497  -0.193   0.356 1.001  2900
## e.y[51]    2.526   0.551   1.438   2.159   2.527   2.893   3.606 1.005   480
## e.y[52]    0.303   0.589  -0.840  -0.100   0.297   0.710   1.452 1.003   760
## e.y[53]   -2.621   0.662  -3.916  -3.074  -2.630  -2.167  -1.315 1.008   290
## e.y[54]    0.459   0.573  -0.637   0.066   0.452   0.851   1.579 1.003  1000
## m1         0.460   0.838  -0.336   0.340   0.461   0.575   1.100 1.081  5200
## m2         2.207   4.589  -2.868   1.521   2.166   2.795   7.467 1.066  2700
## s1         0.389   1.156   0.026   0.078   0.140   0.303   2.484 1.004   770
## s2         3.023   7.570   0.047   0.432   1.043   2.516  19.919 1.001  4200
## sigma.y    1.630   0.171   1.340   1.506   1.615   1.734   2.008 1.001  5200
## deviance 204.454   4.014 198.692 201.449 203.766 206.802 213.816 1.001  5200
## 
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## 
## DIC info (using the rule: pV = var(deviance)/2)
## pV = 8.1 and DIC = 212.5
## DIC is an estimate of expected predictive error (lower deviance is better).
# Calculate R²
SP_model4 <- jags_model_results$BUGSoutput$sims.list   # Extract simulations from the results
r_y_model4 <- 1 - mean(apply(SP_model4$e.y , 1, var)) / var(data$ICA)   # Calculate R² 
cat(paste("R² for model 4:", round(r_y_model4, 4)), "\n")
## R² for model 4: 0.5247
# Extract posterior samples and calculate omega
posterior_samples <- jags_model_results$BUGSoutput$sims.list
omega_model4 <- pmin((apply(posterior_samples$e.y, 2, sd) / mean(posterior_samples$sigma.y))^2, 1)
cat("Omega:\n")
## Omega:
plot(
  omega_model4,
  type = "l", col = "blue", lwd = 2,
  main = "", xlab = "Index", ylab = "Omega"
)

# Define the parameters to analyze
parameters <- c("beta0", "beta1[1]", "beta1[2]", "beta1[3]", "beta2[1]", "beta2[2]", "beta2[3]")

 # Generate and display diagnostic plots 
for (param in parameters) {
  diag_mcmc(as.mcmc(jags_model_results), parName = param)
}